Interior_Mangroove_Sediments
16S Amplicon Microbiota Diversity
Diana Oaxaca
2024-04
Data Prepare
Load data
Check samples distribution
dim(metadata)## [1] 50 30library(ggplot2)
library(ggsci)
# plot
samples_plot <- ggplot(metadata, aes(x = ID, y = ..count.., fill = Depth_cm)) +
geom_bar(position = "dodge") +
facet_wrap(~ Site, scales = "free_x") +
labs(title = "Samples by Site and Depth", y = "Number of Samples", x = "Site") +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) + # Rotate labels
scale_fill_lancet()
#show
samples_plot## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Create Phyloseq
library(phyloseq)
library(qiime2R)
#create phyloseq object
ps <- qza_to_phyloseq(
features = "../results/04.qiime/ASV_table_filter_freq218_emc.qza",
#tree = "../results/04.qiime/rooted-tree-iqtree.qza",
taxonomy = "../results/04.qiime/taxonomy.qza",
metadata = "../data/metadata.tsv")ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]Clean taxonomy table
library(devtools)## Loading required package: usethislibrary(microbiome)##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2022 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>##
## Attaching package: 'microbiome'## The following object is masked from 'package:ggplot2':
##
## alpha## The following object is masked from 'package:base':
##
## transformlibrary(microbiomeutilities)## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'##
## Attaching package: 'microbiomeutilities'## The following object is masked from 'package:microbiome':
##
## add_refseq#check data
microbiome::summarize_phyloseq(ps)## Compositional = NO2## 1] Min. number of reads = 664132] Max. number of reads = 6748463] Total number of reads = 90064984] Average number of reads = 180129.965] Median number of reads = 171486.57] Sparsity = 0.5668845991263916] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 66413"
##
## [[2]]
## [1] "2] Max. number of reads = 674846"
##
## [[3]]
## [1] "3] Total number of reads = 9006498"
##
## [[4]]
## [1] "4] Average number of reads = 180129.96"
##
## [[5]]
## [1] "5] Median number of reads = 171486.5"
##
## [[6]]
## [1] "7] Sparsity = 0.566884599126391"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."Search Phylum with low prevalence
# Explore prevalence
## Get prevalence
prevdf = apply(X = otu_table(ps),
MARGIN = ifelse(taxa_are_rows(ps), yes = 1, no = 2),
FUN = function(x){sum(x > 0)})
## Add taxonomy
prevdf = data.frame(Prevalence = prevdf,
TotalAbundance = taxa_sums(ps),
tax_table(ps))
## Check prevalence at Phylum level
plyr::ddply(prevdf, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})plyr::ddply(prevdf, "Genus", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})## Get genus prevalence plot
prevalence_genus = subset(prevdf, Genus %in% get_taxa_unique(ps, "Genus"))
prev_genus <- ggplot(prevalence_genus, aes(TotalAbundance,
Prevalence /nsamples(ps),color=Genus)) +
geom_hline(yintercept = 0.05, alpha = 0.5, linetype = 2) +
geom_point(size = 1.5, alpha = 0.7) + theme_bw() +
scale_x_log10() + xlab("Total Abundance") + ylab("Prevalence") +
facet_wrap(~Phylum) + theme(legend.position="none")
prev_genus
Clean names
#check names
head(phyloseq::tax_table(ps))## Taxonomy Table: [6 taxa by 7 taxonomic ranks]:
## Kingdom Phylum Class
## ASV_2303 "d__Bacteria" "Proteobacteria" "Alphaproteobacteria"
## ASV_6775 "d__Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_2543 "d__Bacteria" "Chloroflexi" "Dehalococcoidia"
## ASV_4404 "d__Bacteria" "Chloroflexi" "Anaerolineae"
## ASV_1529 "d__Bacteria" "Latescibacterota" "Latescibacterota"
## ASV_779 "d__Bacteria" "Chloroflexi" "Anaerolineae"
## Order Family Genus
## ASV_2303 "Rhizobiales" "Xanthobacteraceae" "Pseudolabrys"
## ASV_6775 "Burkholderiales" "Sutterellaceae" "uncultured"
## ASV_2543 "MSBL5" "MSBL5" "MSBL5"
## ASV_4404 "SBR1031" "SBR1031" "SBR1031"
## ASV_1529 "Latescibacterota" "Latescibacterota" "Latescibacterota"
## ASV_779 "Anaerolineales" "Anaerolineaceae" "uncultured"
## Species
## ASV_2303 NA
## ASV_6775 "uncultured_bacterium"
## ASV_2543 "uncultured_bacterium"
## ASV_4404 "uncultured_bacterium"
## ASV_1529 NA
## ASV_779 "uncultured_Chloroflexi"tail(phyloseq::tax_table(ps))## Taxonomy Table: [6 taxa by 7 taxonomic ranks]:
## Kingdom Phylum Class
## ASV_4505 "d__Bacteria" "MBNT15" "MBNT15"
## ASV_5666 "d__Bacteria" "Chloroflexi" "Dehalococcoidia"
## ASV_108 "d__Bacteria" "Nitrospirota" "Thermodesulfovibrionia"
## ASV_3882 "d__Bacteria" "Firmicutes" "Bacilli"
## ASV_6154 "d__Bacteria" "Methylomirabilota" "Methylomirabilia"
## ASV_1246 "d__Bacteria" "Proteobacteria" "Gammaproteobacteria"
## Order Family Genus
## ASV_4505 "MBNT15" "MBNT15" "MBNT15"
## ASV_5666 "t0.6.f" "t0.6.f" "t0.6.f"
## ASV_108 "uncultured" "uncultured" "uncultured"
## ASV_3882 "Lactobacillales" "Enterococcaceae" "Enterococcus"
## ASV_6154 "Rokubacteriales" "Rokubacteriales" "Rokubacteriales"
## ASV_1246 "Burkholderiales" "Comamonadaceae" "Delftia"
## Species
## ASV_4505 "uncultured_prokaryote"
## ASV_5666 "uncultured_Chloroflexi"
## ASV_108 "uncultured_Nitrospirae"
## ASV_3882 NA
## ASV_6154 "uncultured_bacterium"
## ASV_1246 NAAgglomerate levels
#Aggregate at genus level.
ps_family <- phyloseq::tax_glom(ps, "Family", NArm = TRUE)
ps_genus <- phyloseq::tax_glom(ps, "Genus", NArm = TRUE)
ps_species <- phyloseq::tax_glom(ps, "Species", NArm = TRUE)ps_genus <- filter_taxa(ps_genus, function(x) sum(x) > 1, TRUE)
head(phyloseq::tax_table(ps_genus))## Taxonomy Table: [6 taxa by 7 taxonomic ranks]:
## Kingdom Phylum Class
## ASV_3181 "d__Bacteria" "Myxococcota" "Polyangia"
## ASV_4916 "d__Archaea" "Nanoarchaeota" "Nanoarchaeia"
## ASV_1765 "d__Bacteria" "Planctomycetota" "Phycisphaerae"
## ASV_1959 "d__Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_7171 "d__Bacteria" "Armatimonadota" "Fimbriimonadia"
## ASV_2739 "d__Bacteria" "Verrucomicrobiota" "Verrucomicrobiae"
## Order Family
## ASV_3181 "VHS-B3-70" "VHS-B3-70"
## ASV_4916 "Woesearchaeales" "Woesearchaeales"
## ASV_1765 "Phycisphaerales" "Phycisphaeraceae"
## ASV_1959 "Gammaproteobacteria_Incertae_Sedis" "Unknown_Family"
## ASV_7171 "Fimbriimonadales" "Fimbriimonadaceae"
## ASV_2739 "Opitutales" "Opitutaceae"
## Genus Species
## ASV_3181 "VHS-B3-70" NA
## ASV_4916 "Woesearchaeales" NA
## ASV_1765 "AKYG587" NA
## ASV_1959 "uncultured" NA
## ASV_7171 "Fimbriimonadaceae" NA
## ASV_2739 "Opitutus" NAChange id by tax name
#Substitute these IDs with names of genus
ps_genus <- format_to_besthit(ps_genus) #no me gusta el formato, pero sirvehead(phyloseq::tax_table(ps_genus))## Taxonomy Table: [6 taxa by 8 taxonomic ranks]:
## Domain Phylum
## ASV_3181:g__VHS-B3-70 "d__Bacteria" "Myxococcota"
## ASV_4916:g__Woesearchaeales "d__Archaea" "Nanoarchaeota"
## ASV_1765:g__AKYG587 "d__Bacteria" "Planctomycetota"
## ASV_1959:g__uncultured "d__Bacteria" "Proteobacteria"
## ASV_7171:g__Fimbriimonadaceae "d__Bacteria" "Armatimonadota"
## ASV_2739:g__Opitutus "d__Bacteria" "Verrucomicrobiota"
## Class
## ASV_3181:g__VHS-B3-70 "Polyangia"
## ASV_4916:g__Woesearchaeales "Nanoarchaeia"
## ASV_1765:g__AKYG587 "Phycisphaerae"
## ASV_1959:g__uncultured "Gammaproteobacteria"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadia"
## ASV_2739:g__Opitutus "Verrucomicrobiae"
## Order
## ASV_3181:g__VHS-B3-70 "VHS-B3-70"
## ASV_4916:g__Woesearchaeales "Woesearchaeales"
## ASV_1765:g__AKYG587 "Phycisphaerales"
## ASV_1959:g__uncultured "Gammaproteobacteria_Incertae_Sedis"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadales"
## ASV_2739:g__Opitutus "Opitutales"
## Family Genus
## ASV_3181:g__VHS-B3-70 "VHS-B3-70" "VHS-B3-70"
## ASV_4916:g__Woesearchaeales "Woesearchaeales" "Woesearchaeales"
## ASV_1765:g__AKYG587 "Phycisphaeraceae" "AKYG587"
## ASV_1959:g__uncultured "Unknown_Family" "uncultured"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadaceae" "Fimbriimonadaceae"
## ASV_2739:g__Opitutus "Opitutaceae" "Opitutus"
## Species
## ASV_3181:g__VHS-B3-70 "g__VHS-B3-70"
## ASV_4916:g__Woesearchaeales "g__Woesearchaeales"
## ASV_1765:g__AKYG587 "g__AKYG587"
## ASV_1959:g__uncultured "g__uncultured"
## ASV_7171:g__Fimbriimonadaceae "g__Fimbriimonadaceae"
## ASV_2739:g__Opitutus "g__Opitutus"
## best_hit
## ASV_3181:g__VHS-B3-70 "ASV_3181:g__VHS-B3-70"
## ASV_4916:g__Woesearchaeales "ASV_4916:g__Woesearchaeales"
## ASV_1765:g__AKYG587 "ASV_1765:g__AKYG587"
## ASV_1959:g__uncultured "ASV_1959:g__uncultured"
## ASV_7171:g__Fimbriimonadaceae "ASV_7171:g__Fimbriimonadaceae"
## ASV_2739:g__Opitutus "ASV_2739:g__Opitutus"Clean rare characters
# Remove characters
taxa_names(ps_genus) <- gsub(" ", ".", taxa_names(ps_genus))
taxa_names(ps_genus) <- gsub("\\(|\\)", "", taxa_names(ps_genus))
#check
unique(tax_table(ps_genus)[,"Genus"] )## Taxonomy Table: [462 taxa by 1 taxonomic ranks]:
## Genus
## ASV_3181:g__VHS-B3-70 "VHS-B3-70"
## ASV_4916:g__Woesearchaeales "Woesearchaeales"
## ASV_1765:g__AKYG587 "AKYG587"
## ASV_1959:g__uncultured "uncultured"
## ASV_7171:g__Fimbriimonadaceae "Fimbriimonadaceae"
## ASV_2739:g__Opitutus "Opitutus"
## ASV_2126:g__Cloacibacterium "Cloacibacterium"
## ASV_4050:g__Subgroup_12 "Subgroup_12"
## ASV_3872:g__Hadarchaeales "Hadarchaeales"
## ASV_3878:g__Desulfuromonadaceae "Desulfuromonadaceae"
## ASV_1608:g__DEV008 "DEV008"
## ASV_5186:g__Sulfuritalea "Sulfuritalea"
## ASV_3275:g__Urania-1B-19_marine_sediment_group "Urania-1B-19_marine_sediment_group"
## ASV_840:g__WS2 "WS2"
## ASV_3640:g__Sediminispirochaeta "Sediminispirochaeta"
## ASV_3555:g__Anaerococcus "Anaerococcus"
## ASV_4206:g__Pelomonas "Pelomonas"
## ASV_98:g__Bacillus "Bacillus"
## ASV_702:g__Sulfuricurvum "Sulfuricurvum"
## ASV_1309:g__Comamonas "Comamonas"
## ASV_4053:g__Candidatus_Xiphinematobacter "Candidatus_Xiphinematobacter"
## ASV_5353:g__Amphiplicatus "Amphiplicatus"
## ASV_1850:g__Subgroup_23 "Subgroup_23"
## ASV_541:g__Thiobacillus "Thiobacillus"
## ASV_6812:g__Candidatus_Komeilibacteria "Candidatus_Komeilibacteria"
## ASV_6217:g__RB41 "RB41"
## ASV_2551:g__Curvibacter "Curvibacter"
## ASV_1166:g__Dadabacteriales "Dadabacteriales"
## ASV_511:g__Burkholderia-Caballeronia-Paraburkholderia "Burkholderia-Caballeronia-Paraburkholderia"
## ASV_2769:g__Kazania "Kazania"
## ASV_3171:g__UBA9983 "UBA9983"
## ASV_1838:g__Sphingomonas "Sphingomonas"
## ASV_2758:g__AKYH767 "AKYH767"
## ASV_5014:g__Pseudoclavibacter "Pseudoclavibacter"
## ASV_4921:g__UTBCD1 "UTBCD1"
## ASV_4553:g__[Cytophaga]_xylanolytica_group "[Cytophaga]_xylanolytica_group"
## ASV_3865:g__Kapabacteriales "Kapabacteriales"
## ASV_1920:g__CK-2C2-2 "CK-2C2-2"
## ASV_199:g__KF-JG30-B3 "KF-JG30-B3"
## ASV_1164:g__Chitinibacter "Chitinibacter"
## ASV_183:g__PLTA13 "PLTA13"
## ASV_5281:g__Dyadobacter "Dyadobacter"
## ASV_5:g__Aeromonas "Aeromonas"
## ASV_1384:g__DS-100 "DS-100"
## ASV_89:g__GIF9 "GIF9"
## ASV_3631:g__SCGC_AAA011-D5 "SCGC_AAA011-D5"
## ASV_5107:g__Succinivibrio "Succinivibrio"
## ASV_5917:g__Bacteroidales "Bacteroidales"
## ASV_5408:g__Dinghuibacter "Dinghuibacter"
## ASV_1880:g__Dechlorobacter "Dechlorobacter"
## ASV_1772:g__MB-A2-108 "MB-A2-108"
## ASV_1243:g__Acidibacter "Acidibacter"
## ASV_6811:g__Brevibacillus "Brevibacillus"
## ASV_730:g__LCP-89 "LCP-89"
## ASV_4530:g__Candidatus_Levybacteria "Candidatus_Levybacteria"
## ASV_4582:g__hgcI_clade "hgcI_clade"
## ASV_3225:g__Sphingobium "Sphingobium"
## ASV_3205:g__RBG-16-58-14 "RBG-16-58-14"
## ASV_146:g__Leptospirillum "Leptospirillum"
## ASV_5819:g__A4b "A4b"
## ASV_501:g__Proteus "Proteus"
## ASV_1065:g__SCGC_AB-179-E04 "SCGC_AB-179-E04"
## ASV_1288:g__Rhodoferax "Rhodoferax"
## ASV_1460:g__Anaeromyxobacter "Anaeromyxobacter"
## ASV_260:g__Subgroup_9 "Subgroup_9"
## ASV_1313:g__ADurb.Bin180 "ADurb.Bin180"
## ASV_997:g__Vermiphilaceae "Vermiphilaceae"
## ASV_6192:g__Lineage_IV "Lineage_IV"
## ASV_1909:g__KD4-96 "KD4-96"
## ASV_401:g__Calditrichaceae "Calditrichaceae"
## ASV_113:g__Zixibacteria "Zixibacteria"
## ASV_670:g__SC-I-84 "SC-I-84"
## ASV_495:g__RBG-13-46-9 "RBG-13-46-9"
## ASV_2598:g__10bav-F6 "10bav-F6"
## ASV_2545:g__37-13 "37-13"
## ASV_201:g__Entomoplasma "Entomoplasma"
## ASV_235:g__ADurb.Bin118 "ADurb.Bin118"
## ASV_470:g__Weissella "Weissella"
## ASV_393:g__TSAC18 "TSAC18"
## ASV_2272:g__S0134_terrestrial_group "S0134_terrestrial_group"
## ASV_84:g__Gemmobacter "Gemmobacter"
## ASV_502:g__Nitrospira "Nitrospira"
## ASV_139:g__661239 "661239"
## ASV_5789:g__Ochrobactrum "Ochrobactrum"
## ASV_383:g__Asteroleplasma "Asteroleplasma"
## ASV_2215:g__Pantoea "Pantoea"
## ASV_1264:g__Leptotrichia "Leptotrichia"
## ASV_4:g__Bifidobacterium "Bifidobacterium"
## ASV_5453:g__Chthoniobacter "Chthoniobacter"
## ASV_2401:g__Fictibacillus "Fictibacillus"
## ASV_2581:g__WOR-1 "WOR-1"
## ASV_122:g__Rheinheimera "Rheinheimera"
## ASV_2339:g__Aquabacterium "Aquabacterium"
## ASV_266:g__Hirschia "Hirschia"
## ASV_227:g__Staphylococcus "Staphylococcus"
## ASV_171:g__MSBL5 "MSBL5"
## ASV_5898:g__Agathobacter "Agathobacter"
## ASV_446:g__Clostridium_sensu_stricto_8 "Clostridium_sensu_stricto_8"
## ASV_1757:g__OLB14 "OLB14"
## ASV_4226:g__PS-B29 "PS-B29"
## ASV_97:g__Morganella "Morganella"
## ASV_5731:g__Obscuribacteraceae "Obscuribacteraceae"
## ASV_296:g__ER-E4-17 "ER-E4-17"
## ASV_2238:g__Chryseolinea "Chryseolinea"
## ASV_587:g__Subgroup_2 "Subgroup_2"
## ASV_5651:g__Sulfuriferula "Sulfuriferula"
## ASV_886:g__Dechloromonas "Dechloromonas"
## ASV_2301:g__Streptomyces "Streptomyces"
## ASV_1093:g__Providencia "Providencia"
## ASV_1377:g__Candidatus_Collierbacteria "Candidatus_Collierbacteria"
## ASV_1646:g__UCG-004 "UCG-004"
## ASV_596:g__Subgroup_17 "Subgroup_17"
## ASV_352:g__mle1-7 "mle1-7"
## ASV_134:g__Desulfobacca "Desulfobacca"
## ASV_1992:g__AKAU3564_sediment_group "AKAU3564_sediment_group"
## ASV_1091:g__MBMPE27 "MBMPE27"
## ASV_1785:g__Pajaroellobacter "Pajaroellobacter"
## ASV_463:g__SJA-28 "SJA-28"
## ASV_111:g__MBNT15 "MBNT15"
## ASV_276:g__BSV26 "BSV26"
## ASV_6374:g__Desulfatirhabdium "Desulfatirhabdium"
## ASV_59:g__Pectinatus "Pectinatus"
## ASV_563:g__GOUTB8 "GOUTB8"
## ASV_5086:g__Acidobacteriota "Acidobacteriota"
## ASV_4953:g__Thermodesulfovibrio "Thermodesulfovibrio"
## ASV_1516:g__Kineosporia "Kineosporia"
## ASV_58:g__Bradyrhizobium "Bradyrhizobium"
## ASV_683:g__Syntrophorhabdus "Syntrophorhabdus"
## ASV_1030:g__bacteriap25 "bacteriap25"
## ASV_414:g__SAR202_clade "SAR202_clade"
## ASV_7033:g__Micromonospora "Micromonospora"
## ASV_2640:g__SJA-15 "SJA-15"
## ASV_3609:g__Paucibacter "Paucibacter"
## ASV_3740:g__Brevundimonas "Brevundimonas"
## ASV_3379:g__Lawsonella "Lawsonella"
## ASV_2528:g__Marinimicrobia_SAR406_clade "Marinimicrobia_(SAR406_clade)"
## ASV_516:g__Nocardiopsis "Nocardiopsis"
## ASV_4827:g__Candidatus_Doudnabacteria "Candidatus_Doudnabacteria"
## ASV_126:g__H3.93 "H3.93"
## ASV_5109:g__Gemmata "Gemmata"
## ASV_1401:g__PAUC26f "PAUC26f"
## ASV_1936:g__Cupriavidus "Cupriavidus"
## ASV_4746:g__Treponema "Treponema"
## ASV_644:g__GAL15 "GAL15"
## ASV_1967:g__Candidatus_Kaiserbacteria "Candidatus_Kaiserbacteria"
## ASV_4485:g__Vagococcus "Vagococcus"
## ASV_204:g__Methylobacterium-Methylorubrum "Methylobacterium-Methylorubrum"
## ASV_3596:g__Candidatus_Moranbacteria "Candidatus_Moranbacteria"
## ASV_1903:g__Z195MB87 "Z195MB87"
## ASV_350:g__Clostridium_sensu_stricto_5 "Clostridium_sensu_stricto_5"
## ASV_3246:g__IMCC26256 "IMCC26256"
## ASV_877:g__Flavobacterium "Flavobacterium"
## ASV_67:g__Spirochaeta "Spirochaeta"
## ASV_92:g__GIF3 "GIF3"
## ASV_465:g__Sh765B-TzT-20 "Sh765B-TzT-20"
## ASV_2317:g__Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"
## ASV_6839:g__P3OB-42 "P3OB-42"
## ASV_599:g__Chitinimonas "Chitinimonas"
## ASV_1092:g__Marine_Benthic_Group_D_and_DHVEG-1 "Marine_Benthic_Group_D_and_DHVEG-1"
## ASV_1392:g__Romboutsia "Romboutsia"
## ASV_1256:g__Sericytochromatia "Sericytochromatia"
## ASV_3182:g__Ga0077536 "Ga0077536"
## ASV_148:g__S085 "S085"
## ASV_226:g__Streptococcus "Streptococcus"
## ASV_429:g__Hypnocyclicus "Hypnocyclicus"
## ASV_356:g__GOUTA6 "GOUTA6"
## ASV_6978:g__Hafnia-Obesumbacterium "Hafnia-Obesumbacterium"
## ASV_6342:g__Gracilibacteria "Gracilibacteria"
## ASV_32:g__Aminicenantales "Aminicenantales"
## ASV_1238:g__SPG12-343-353-B75 "SPG12-343-353-B75"
## ASV_2547:g__Lacunisphaera "Lacunisphaera"
## ASV_769:g__PHOS-HE36 "PHOS-HE36"
## ASV_4642:g__MSB-5E12 "MSB-5E12"
## ASV_5286:g__Dyella "Dyella"
## ASV_362:g__Dysgonomonas "Dysgonomonas"
## ASV_1519:g__Subgroup_22 "Subgroup_22"
## ASV_132:g__SB-5 "SB-5"
## ASV_6145:g__Pedosphaera "Pedosphaera"
## ASV_1140:g__AKIW659 "AKIW659"
## ASV_925:g__SEEP-SRB1 "SEEP-SRB1"
## ASV_2297:g__Schlegelella "Schlegelella"
## ASV_5112:g__Fimbriiglobus "Fimbriiglobus"
## ASV_6393:g__Phyllobacterium "Phyllobacterium"
## ASV_852:g__FS117-23B-02 "FS117-23B-02"
## ASV_3317:g__Gitt-GS-136 "Gitt-GS-136"
## ASV_2927:g__Hydrogenedensaceae "Hydrogenedensaceae"
## ASV_5128:g__Dehalogenimonas "Dehalogenimonas"
## ASV_1637:g__LD1-PA32 "LD1-PA32"
## ASV_1492:g__CCM11a "CCM11a"
## ASV_1220:g__CG2-30-50-142 "CG2-30-50-142"
## ASV_6132:g__Fusobacterium "Fusobacterium"
## ASV_160:g__KF-JG30-C25 "KF-JG30-C25"
## ASV_1308:g__S-BQ2-57_soil_group "S-BQ2-57_soil_group"
## ASV_3908:g__LWQ8 "LWQ8"
## ASV_223:g__BD2-11_terrestrial_group "BD2-11_terrestrial_group"
## ASV_2414:g__Deferrisoma "Deferrisoma"
## ASV_21:g__Micrococcus "Micrococcus"
## ASV_4595:g__Clostridium_sensu_stricto_9 "Clostridium_sensu_stricto_9"
## ASV_585:g__Subgroup_7 "Subgroup_7"
## ASV_1851:g__BSN166 "BSN166"
## ASV_1702:g__WD2101_soil_group "WD2101_soil_group"
## ASV_297:g__Calorithrix "Calorithrix"
## ASV_5723:g__Pir1_lineage "Pir1_lineage"
## ASV_4999:g__Run-SP154 "Run-SP154"
## ASV_1564:g__Thioalkalispira-Sulfurivermis "Thioalkalispira-Sulfurivermis"
## ASV_4344:g__OC31 "OC31"
## ASV_317:g__Megasphaera "Megasphaera"
## ASV_6814:g__CI75cm.2.12 "CI75cm.2.12"
## ASV_27:g__TRA3-20 "TRA3-20"
## ASV_921:g__Oceanobacillus "Oceanobacillus"
## ASV_1124:g__Faecalibacterium "Faecalibacterium"
## ASV_1990:g__Candidatus_Adlerbacteria "Candidatus_Adlerbacteria"
## ASV_834:g__Thermoanaerobaculum "Thermoanaerobaculum"
## ASV_5403:g__Sphingobacterium "Sphingobacterium"
## ASV_320:g__EPR3968-O8a-Bc78 "EPR3968-O8a-Bc78"
## ASV_1009:g__Moduliflexaceae "Moduliflexaceae"
## ASV_1932:g__Candidatus_Peregrinibacteria "Candidatus_Peregrinibacteria"
## ASV_1075:g__Rhodococcus "Rhodococcus"
## ASV_1491:g__mle1-8 "mle1-8"
## ASV_167:g__Deinococcus "Deinococcus"
## ASV_4158:g__WN-HWB-116 "WN-HWB-116"
## ASV_1873:g__Stenotrophomonas "Stenotrophomonas"
## ASV_4902:g__JdFR-54 "JdFR-54"
## ASV_802:g__Raoultella "Raoultella"
## ASV_1249:g__Crenobacter "Crenobacter"
## ASV_1520:g__WX65 "WX65"
## ASV_131:g__A0839 "A0839"
## ASV_221:g__vadinBA26 "vadinBA26"
## ASV_1950:g__Acetothermiia "Acetothermiia"
## ASV_1518:g__Shimwellia "Shimwellia"
## ASV_693:g__RCP2-54 "RCP2-54"
## ASV_4558:g__Luteolibacter "Luteolibacter"
## ASV_5802:g__Desulfovibrio "Desulfovibrio"
## ASV_2653:g__Pediococcus "Pediococcus"
## ASV_325:g__WWE3 "WWE3"
## ASV_3582:g__Schekmanbacteria "Schekmanbacteria"
## ASV_1341:g__Maricaulis "Maricaulis"
## ASV_2727:g__ADurb.Bin063-1 "ADurb.Bin063-1"
## ASV_1779:g__Caldithrix "Caldithrix"
## ASV_194:g__FW22 "FW22"
## ASV_1254:g__Arcobacter "Arcobacter"
## ASV_1684:g__Fimbriimonadales "Fimbriimonadales"
## ASV_6678:g__Lacibacter "Lacibacter"
## ASV_3700:g__B103G10 "B103G10"
## ASV_1919:g__Pir4_lineage "Pir4_lineage"
## ASV_2957:g__Pedosphaeraceae "Pedosphaeraceae"
## ASV_5692:g__Nordella "Nordella"
## ASV_286:g__Lactococcus "Lactococcus"
## ASV_5697:g__Steroidobacter "Steroidobacter"
## ASV_3155:g__Z-35 "Z-35"
## ASV_4374:g__Moheibacter "Moheibacter"
## ASV_1691:g__Desulfomonile "Desulfomonile"
## ASV_2342:g__Syntrophobacter "Syntrophobacter"
## ASV_26:g__Napoli-4B-65 "Napoli-4B-65"
## ASV_2928:g__Sandaracinus "Sandaracinus"
## ASV_1818:g__Domibacillus "Domibacillus"
## ASV_705:g__Bauldia "Bauldia"
## ASV_24:g__Escherichia-Shigella "Escherichia-Shigella"
## ASV_2434:g__Pla4_lineage "Pla4_lineage"
## ASV_5301:g__Megamonas "Megamonas"
## ASV_4731:g__Actinomyces "Actinomyces"
## ASV_118:g__MND1 "MND1"
## ASV_1862:g__Subgroup_25 "Subgroup_25"
## ASV_2521:g__Candidatus_Liptonbacteria "Candidatus_Liptonbacteria"
## ASV_6258:g__WS4 "WS4"
## ASV_477:g__Subgroup_15 "Subgroup_15"
## ASV_1615:g__B2M28 "B2M28"
## ASV_2747:g__Cutibacterium "Cutibacterium"
## ASV_6917:g__Chitinophaga "Chitinophaga"
## ASV_4295:g__Subgroup_10 "Subgroup_10"
## ASV_3212:g__BIrii41 "BIrii41"
## ASV_1276:g__Gammaproteobacteria "Gammaproteobacteria"
## ASV_208:g__Dehalococcoidia "Dehalococcoidia"
## ASV_5788:g__Candidatus_Uhrbacteria "Candidatus_Uhrbacteria"
## ASV_7160:g__Pla3_lineage "Pla3_lineage"
## ASV_3495:g__Blfdi19 "Blfdi19"
## ASV_1489:g__Candidatus_Magasanikbacteria "Candidatus_Magasanikbacteria"
## ASV_5104:g__Thiogranum "Thiogranum"
## ASV_1027:g__ODP1230B23.02 "ODP1230B23.02"
## ASV_4874:g__Eel-36e1D6 "Eel-36e1D6"
## ASV_115:g__Bathyarchaeia "Bathyarchaeia"
## ASV_1325:g__Terrimonas "Terrimonas"
## ASV_2644:g__Clostridium_sensu_stricto_13 "Clostridium_sensu_stricto_13"
## ASV_105:g__Sva0081_sediment_group "Sva0081_sediment_group"
## ASV_7:g__Sva0485 "Sva0485"
## ASV_137:g__Exiguobacterium "Exiguobacterium"
## ASV_374:g__Thermoflexus "Thermoflexus"
## ASV_3438:g__oc32 "oc32"
## ASV_1517:g__JG36-GS-52 "JG36-GS-52"
## ASV_3269:g__Gemmatimonadaceae "Gemmatimonadaceae"
## ASV_678:g__RBG-16-49-21 "RBG-16-49-21"
## ASV_5555:g__Shewanella "Shewanella"
## ASV_6776:g__Bradymonadales "Bradymonadales"
## ASV_22:g__Lactobacillus "Lactobacillus"
## ASV_3261:g__SH-PL14 "SH-PL14"
## ASV_192:g__4-29-1 "4-29-1"
## ASV_709:g__Subgroup_13 "Subgroup_13"
## ASV_3165:g__Rhodomicrobium "Rhodomicrobium"
## ASV_738:g__Bacteroidetes_BD2-2 "Bacteroidetes_BD2-2"
## ASV_456:g__SBR1031 "SBR1031"
## ASV_1642:g__Candidatus_Paceibacter "Candidatus_Paceibacter"
## ASV_6714:g__Proteiniphilum "Proteiniphilum"
## ASV_593:g__Sporacetigenium "Sporacetigenium"
## ASV_198:g__Serratia "Serratia"
## ASV_1797:g__Odinarchaeia "Odinarchaeia"
## ASV_4836:g__Candidatus_Alysiosphaera "Candidatus_Alysiosphaera"
## ASV_3711:g__Chitinilyticum "Chitinilyticum"
## ASV_6655:g__Subdoligranulum "Subdoligranulum"
## ASV_3402:g__Ruminiclostridium "Ruminiclostridium"
## ASV_3089:g__Altererythrobacter "Altererythrobacter"
## ASV_4278:g__LD-RB-34 "LD-RB-34"
## ASV_200:g__MD2902-B12 "MD2902-B12"
## ASV_6053:g__Holophagaceae "Holophagaceae"
## ASV_6512:g__UCG-012 "UCG-012"
## ASV_5879:g__Ferruginibacter "Ferruginibacter"
## ASV_2070:g__SG8-4 "SG8-4"
## ASV_55:g__Acinetobacter "Acinetobacter"
## ASV_272:g__Pedomicrobium "Pedomicrobium"
## ASV_51:g__Enhydrobacter "Enhydrobacter"
## ASV_459:g__Subgroup_21 "Subgroup_21"
## ASV_3096:g__KI89A_clade "KI89A_clade"
## ASV_5622:g__Fervidicella "Fervidicella"
## ASV_4617:g__Candidatus_Chisholmbacteria "Candidatus_Chisholmbacteria"
## ASV_265:g__Latescibacteraceae "Latescibacteraceae"
## ASV_1011:g__DG-56 "DG-56"
## ASV_110:g__ZOR0006 "ZOR0006"
## ASV_4811:g__Candidatus_Thiosymbion "Candidatus_Thiosymbion"
## ASV_3791:g__A21b "A21b"
## ASV_1189:g__Luteimonas "Luteimonas"
## ASV_240:g__Sh765B-AG-111 "Sh765B-AG-111"
## ASV_990:g__Dongia "Dongia"
## ASV_6472:g__Rivicola "Rivicola"
## ASV_2371:g__SHA-26 "SHA-26"
## ASV_2534:g__Tumebacillus "Tumebacillus"
## ASV_1099:g__Hyphomicrobium "Hyphomicrobium"
## ASV_2232:g__Inhella "Inhella"
## ASV_3399:g__P9X2b3D02 "P9X2b3D02"
## ASV_2599:g__Candidatus_Latescibacter "Candidatus_Latescibacter"
## ASV_2105:g__Candidatus_Solibacter "Candidatus_Solibacter"
## ASV_4958:g__Woeseia "Woeseia"
## ASV_5875:g__Entotheonellaceae "Entotheonellaceae"
## ASV_1465:g__Clostridium_sensu_stricto_6 "Clostridium_sensu_stricto_6"
## ASV_1991:g__Subgroup_5 "Subgroup_5"
## ASV_4177:g__Limibaculum "Limibaculum"
## ASV_5483:g__[Polyangium]_brachysporum_group "[Polyangium]_brachysporum_group"
## ASV_73:g__t0.6.f "t0.6.f"
## ASV_69:g__Desulfatiglans "Desulfatiglans"
## ASV_2235:g__Ensifer "Ensifer"
## ASV_5272:g__Candidatus_Portnoybacteria "Candidatus_Portnoybacteria"
## ASV_5890:g__Candidatus_Methylospira "Candidatus_Methylospira"
## ASV_4434:g__Nitrolancea "Nitrolancea"
## ASV_4377:g__Pla1_lineage "Pla1_lineage"
## ASV_360:g__NB1-j "NB1-j"
## ASV_70:g__SCGC-AB-539-J10 "SCGC-AB-539-J10"
## ASV_6034:g__SAR324_cladeMarine_group_B "SAR324_clade(Marine_group_B)"
## ASV_1291:g__Parcubacteria "Parcubacteria"
## ASV_1461:g__vadinHA49 "vadinHA49"
## ASV_2782:g__Saccharimonadales "Saccharimonadales"
## ASV_433:g__Latescibacterota "Latescibacterota"
## ASV_1706:g__AKAU4049 "AKAU4049"
## ASV_4985:g__Novosphingobium "Novosphingobium"
## ASV_5496:g__EC3 "EC3"
## ASV_3068:g__WCHB1-81 "WCHB1-81"
## ASV_141:g__Corynebacterium "Corynebacterium"
## ASV_676:g__Clostridium_sensu_stricto_1 "Clostridium_sensu_stricto_1"
## ASV_2785:g__Haliangium "Haliangium"
## ASV_4788:g__Pirellula "Pirellula"
## ASV_4272:g__Puia "Puia"
## ASV_3:g__Enterobacter "Enterobacter"
## ASV_480:g__TM7a "TM7a"
## ASV_5877:g__Solitalea "Solitalea"
## ASV_313:g__Ellin516 "Ellin516"
## ASV_3257:g__Candidatus_Woesebacteria "Candidatus_Woesebacteria"
## ASV_39:g__Sulfurifustis "Sulfurifustis"
## ASV_450:g__Brachybacterium "Brachybacterium"
## ASV_5386:g__Permianibacter "Permianibacter"
## ASV_4805:g__3PJM14 "3PJM14"
## ASV_5535:g__Leeia "Leeia"
## ASV_519:g__Poribacteria "Poribacteria"
## ASV_3723:g__Amsterdam-1B-07 "Amsterdam-1B-07"
## ASV_34:g__Subgroup_18 "Subgroup_18"
## ASV_2077:g__Lysobacter "Lysobacter"
## ASV_3408:g__wb1-A12 "wb1-A12"
## ASV_228:g__Ignavibacterium "Ignavibacterium"
## ASV_539:g__Pelolinea "Pelolinea"
## ASV_4701:g__UASB-TL25 "UASB-TL25"
## ASV_5499:g__S-70 "S-70"
## ASV_7088:g__EF100-94H03 "EF100-94H03"
## ASV_6730:g__Candidatus_Nomurabacteria "Candidatus_Nomurabacteria"
## ASV_5459:g__MM2 "MM2"
## ASV_1287:g__BD1-7_clade "BD1-7_clade"
## ASV_1:g__Enterococcus "Enterococcus"
## ASV_4188:g__WPS-2 "WPS-2"
## ASV_4223:g__Parapedobacter "Parapedobacter"
## ASV_112:g__Rokubacteriales "Rokubacteriales"
## ASV_857:g__Candidatus_Yanofskybacteria "Candidatus_Yanofskybacteria"
## ASV_5663:g__966-1 "966-1"
## ASV_775:g__Bacteroidetes_vadinHA17 "Bacteroidetes_vadinHA17"
## ASV_1356:g__Ralstonia "Ralstonia"
## ASV_819:g__RBG-13-54-9 "RBG-13-54-9"
## ASV_3952:g__Subgroup_11 "Subgroup_11"
## ASV_746:g__Massilia "Massilia"
## ASV_4792:g__Verrucomicrobium "Verrucomicrobium"
## ASV_6367:g__Candidatus_Jorgensenbacteria "Candidatus_Jorgensenbacteria"
## ASV_3672:g__OM190 "OM190"
## ASV_130:g__Acetobacter "Acetobacter"
## ASV_4124:g__Sideroxydans "Sideroxydans"
## ASV_4111:g__Tolumonas "Tolumonas"
## ASV_752:g__Bryobacter "Bryobacter"
## ASV_703:g__Candidatus_Omnitrophus "Candidatus_Omnitrophus"
## ASV_3320:g__Rhodobacter "Rhodobacter"
## ASV_5092:g__Longimicrobiaceae "Longimicrobiaceae"
## ASV_3517:g__OPB41 "OPB41"
## ASV_66:g__Sh765B-TzT-35 "Sh765B-TzT-35"
## ASV_4666:g__CL500-29_marine_group "CL500-29_marine_group"
## ASV_72:g__Vogesella "Vogesella"
## ASV_3529:g__Thermincola "Thermincola"
## ASV_608:g__CCD24 "CCD24"
## ASV_590:g__Thiobacter "Thiobacter"
## ASV_5373:g__Microgenomatia "Microgenomatia"
## ASV_5790:g__Bacteroidetes_VC2.1_Bac22 "Bacteroidetes_VC2.1_Bac22"
## ASV_41:g__B1-7BS "B1-7BS"
## ASV_942:g__Olsenella "Olsenella"
## ASV_6267:g__Devosia "Devosia"
## ASV_725:g__JG30-KF-CM66 "JG30-KF-CM66"
## ASV_513:g__Sumerlaea "Sumerlaea"
## ASV_2511:g__Reyranella "Reyranella"
## ASV_4756:g__Blastococcus "Blastococcus"
## ASV_1041:g__Mesorhizobium "Mesorhizobium"
## ASV_5164:g__RBG-16-55-12 "RBG-16-55-12"
## ASV_1001:g__AB-539-J10 "AB-539-J10"
## ASV_555:g__Hydrogenophaga "Hydrogenophaga"
## ASV_856:g__Pseudolabrys "Pseudolabrys"
## ASV_184:g__Vicinamibacteraceae "Vicinamibacteraceae"
## ASV_2719:g__SM23-31 "SM23-31"
## ASV_4463:g__OM27_clade "OM27_clade"
## ASV_797:g__Nitratiruptor "Nitratiruptor"
## ASV_1078:g__GWE2-31-10 "GWE2-31-10"
## ASV_3338:g__Paenarthrobacter "Paenarthrobacter"
## ASV_6583:g__DEV114 "DEV114"
## ASV_8:g__Prevotella "Prevotella"
## ASV_3311:g__Nitrosopumilaceae "Nitrosopumilaceae"
## ASV_2322:g__Methyloligellaceae "Methyloligellaceae"
## ASV_1839:g__Leuconostoc "Leuconostoc"
## ASV_5937:g__Azohydromonas "Azohydromonas"
## ASV_1423:g__Gaiella "Gaiella"
## ASV_2282:g__Vulgatibacter "Vulgatibacter"
## ASV_1114:g__Ellin6067 "Ellin6067"
## ASV_1503:g__MSB-5B2 "MSB-5B2"
## ASV_3486:g__Lentimicrobiaceae "Lentimicrobiaceae"
## ASV_812:g__Z114MB74 "Z114MB74"
## ASV_979:g__Edwardsbacteria "Edwardsbacteria"
## ASV_611:g__11-24 "11-24"
## ASV_7101:g__JG36-TzT-191 "JG36-TzT-191"
## ASV_6:g__Pseudomonas "Pseudomonas"
## ASV_935:g__Dialister "Dialister"
## ASV_2933:g__Pseudochrobactrum "Pseudochrobactrum"
## ASV_328:g__TA06 "TA06"
## ASV_355:g__SWB02 "SWB02"
## ASV_1639:g__Fluviicola "Fluviicola"
## ASV_2616:g__Kytococcus "Kytococcus"
## ASV_1246:g__Delftia "Delftia"Subset Sites phyloseq objects
# check data
table(meta(ps)$Site, meta(ps)$Depth_cm)##
## 0-15 16-30 31-45 50-75
## Site_1 2 2 2 2
## Site_2 2 2 2 2
## Site_3 2 2 2 2
## Site_4 1 1 1 1
## Site_5 2 2 2 2
## Site_6 2 2 2 0
## Site_7 2 2 2 2#subsets
ps_s1 <- subset_samples(ps, Site == "Site_1")
ps_s2 <- subset_samples(ps, Site == "Site_2")
ps_s3 <- subset_samples(ps, Site == "Site_3")
ps_s4 <- subset_samples(ps, Site == "Site_4")
ps_s5 <- subset_samples(ps, Site == "Site_5")
ps_s6 <- subset_samples(ps, Site == "Site_6")
ps_s7 <- subset_samples(ps, Site == "Site_7")
#show
ps_s1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]ps_s2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]ps_s3## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]ps_s4## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 4 samples ]
## sample_data() Sample Data: [ 4 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]ps_s5## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]ps_s6## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 6 samples ]
## sample_data() Sample Data: [ 6 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]ps_s7## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]#check data
microbiome::summarize_phyloseq(ps_s1)## Compositional = NO2## 1] Min. number of reads = 1486352] Max. number of reads = 2188773] Total number of reads = 14291624] Average number of reads = 178645.255] Median number of reads = 1712277] Sparsity = 0.5537374947160776] Any OTU sum to 1 or less? YES8] Number of singletons = 7879] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 148635"
##
## [[2]]
## [1] "2] Max. number of reads = 218877"
##
## [[3]]
## [1] "3] Total number of reads = 1429162"
##
## [[4]]
## [1] "4] Average number of reads = 178645.25"
##
## [[5]]
## [1] "5] Median number of reads = 171227"
##
## [[6]]
## [1] "7] Sparsity = 0.553737494716077"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 787"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."microbiome::summarize_phyloseq(ps_s2)## Compositional = NO2## 1] Min. number of reads = 1323482] Max. number of reads = 2205473] Total number of reads = 13798474] Average number of reads = 172480.8755] Median number of reads = 173469.57] Sparsity = 0.5038572636325216] Any OTU sum to 1 or less? YES8] Number of singletons = 12669] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 132348"
##
## [[2]]
## [1] "2] Max. number of reads = 220547"
##
## [[3]]
## [1] "3] Total number of reads = 1379847"
##
## [[4]]
## [1] "4] Average number of reads = 172480.875"
##
## [[5]]
## [1] "5] Median number of reads = 173469.5"
##
## [[6]]
## [1] "7] Sparsity = 0.503857263632521"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 1266"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."microbiome::summarize_phyloseq(ps_s3)## Compositional = NO2## 1] Min. number of reads = 1019792] Max. number of reads = 6748463] Total number of reads = 16868954] Average number of reads = 210861.8755] Median number of reads = 1415347] Sparsity = 0.6287691982527836] Any OTU sum to 1 or less? YES8] Number of singletons = 23169] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 101979"
##
## [[2]]
## [1] "2] Max. number of reads = 674846"
##
## [[3]]
## [1] "3] Total number of reads = 1686895"
##
## [[4]]
## [1] "4] Average number of reads = 210861.875"
##
## [[5]]
## [1] "5] Median number of reads = 141534"
##
## [[6]]
## [1] "7] Sparsity = 0.628769198252783"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 2316"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."microbiome::summarize_phyloseq(ps_s4)## Compositional = NO2## 1] Min. number of reads = 1276312] Max. number of reads = 2065043] Total number of reads = 7156094] Average number of reads = 178902.255] Median number of reads = 1907377] Sparsity = 0.4716077215724956] Any OTU sum to 1 or less? YES8] Number of singletons = 16899] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 127631"
##
## [[2]]
## [1] "2] Max. number of reads = 206504"
##
## [[3]]
## [1] "3] Total number of reads = 715609"
##
## [[4]]
## [1] "4] Average number of reads = 178902.25"
##
## [[5]]
## [1] "5] Median number of reads = 190737"
##
## [[6]]
## [1] "7] Sparsity = 0.471607721572495"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 1689"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."microbiome::summarize_phyloseq(ps_s5)## Compositional = NO2## 1] Min. number of reads = 970812] Max. number of reads = 2982253] Total number of reads = 14629924] Average number of reads = 1828745] Median number of reads = 1679587] Sparsity = 0.5306643652247436] Any OTU sum to 1 or less? YES8] Number of singletons = 8929] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 97081"
##
## [[2]]
## [1] "2] Max. number of reads = 298225"
##
## [[3]]
## [1] "3] Total number of reads = 1462992"
##
## [[4]]
## [1] "4] Average number of reads = 182874"
##
## [[5]]
## [1] "5] Median number of reads = 167958"
##
## [[6]]
## [1] "7] Sparsity = 0.530664365224743"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 892"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."microbiome::summarize_phyloseq(ps_s6)## Compositional = NO2## 1] Min. number of reads = 1549192] Max. number of reads = 2019623] Total number of reads = 10440904] Average number of reads = 1740155] Median number of reads = 172266.57] Sparsity = 0.656803344136026] Any OTU sum to 1 or less? YES8] Number of singletons = 29849] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 154919"
##
## [[2]]
## [1] "2] Max. number of reads = 201962"
##
## [[3]]
## [1] "3] Total number of reads = 1044090"
##
## [[4]]
## [1] "4] Average number of reads = 174015"
##
## [[5]]
## [1] "5] Median number of reads = 172266.5"
##
## [[6]]
## [1] "7] Sparsity = 0.65680334413602"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 2984"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."microbiome::summarize_phyloseq(ps_s7)## Compositional = NO2## 1] Min. number of reads = 664132] Max. number of reads = 2157403] Total number of reads = 12879034] Average number of reads = 160987.8755] Median number of reads = 174806.57] Sparsity = 0.597594053825566] Any OTU sum to 1 or less? YES8] Number of singletons = 22379] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 66413"
##
## [[2]]
## [1] "2] Max. number of reads = 215740"
##
## [[3]]
## [1] "3] Total number of reads = 1287903"
##
## [[4]]
## [1] "4] Average number of reads = 160987.875"
##
## [[5]]
## [1] "5] Median number of reads = 174806.5"
##
## [[6]]
## [1] "7] Sparsity = 0.59759405382556"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
##
## [[8]]
## [1] "8] Number of singletons = 2237"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."Remove singletons
#Delete singletons
ps_s1 <- filter_taxa(ps_s1, function(x) sum(x) > 1, TRUE)
ps_s2 <- filter_taxa(ps_s2, function(x) sum(x) > 1, TRUE)
ps_s3 <- filter_taxa(ps_s3, function(x) sum(x) > 1, TRUE)
ps_s4 <- filter_taxa(ps_s4, function(x) sum(x) > 1, TRUE)
ps_s5 <- filter_taxa(ps_s5, function(x) sum(x) > 1, TRUE)
ps_s6 <- filter_taxa(ps_s6, function(x) sum(x) > 1, TRUE)
ps_s7 <- filter_taxa(ps_s7, function(x) sum(x) > 1, TRUE)
#check
ps_s1## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6310 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6310 taxa by 7 taxonomic ranks ]ps_s2## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5831 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5831 taxa by 7 taxonomic ranks ]ps_s3## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4781 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4781 taxa by 7 taxonomic ranks ]ps_s4## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5408 taxa and 4 samples ]
## sample_data() Sample Data: [ 4 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5408 taxa by 7 taxonomic ranks ]ps_s5## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6205 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6205 taxa by 7 taxonomic ranks ]ps_s6## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4113 taxa and 6 samples ]
## sample_data() Sample Data: [ 6 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4113 taxa by 7 taxonomic ranks ]ps_s7## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4860 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4860 taxa by 7 taxonomic ranks ]Relative abundance phyloseq objects
library(microbiome)
library(microbiomeutilities)
#function
convert_to_compositional_and_index_reformat <- function(ps) {
ps_rel <- microbiome::transform(ps, "compositional")
ps_rel_f <- format_to_besthit(ps_rel)
return(ps_rel_f)
}
# phyloseq objects
phyloseq_objects <- list(ps_s1 = ps_s1, ps_s2 = ps_s2, ps_s3 = ps_s3,
ps_s4 = ps_s4, ps_s5 = ps_s5, ps_s6 = ps_s6,
ps_s7 = ps_s7, ps=ps)
# Apply
phyloseq_objects_rel <- lapply(phyloseq_objects, convert_to_compositional_and_index_reformat)
# Subset
ps_s1_rel_f <- phyloseq_objects_rel$ps_s1
ps_s2_rel_f <- phyloseq_objects_rel$ps_s2
ps_s3_rel_f <- phyloseq_objects_rel$ps_s3
ps_s4_rel_f <- phyloseq_objects_rel$ps_s4
ps_s5_rel_f <- phyloseq_objects_rel$ps_s5
ps_s6_rel_f <- phyloseq_objects_rel$ps_s6
ps_s7_rel_f <- phyloseq_objects_rel$ps_s7
ps_full_rel_f <- phyloseq_objects_rel$psFull core object
# Get full core with microbiome
full_ps_core_relative <- core(ps_full_rel_f, 0, 0.5)
full_ps_core_counts <- core(ps, 0, 0.5)
# Check point
ps## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]ps_full_rel_f## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 8 taxonomic ranks ]full_ps_core_relative## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2048 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 2048 taxa by 8 taxonomic ranks ]full_ps_core_counts## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2048 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 2048 taxa by 7 taxonomic ranks ]microbiome::summarize_phyloseq(full_ps_core_counts)## Compositional = NO2## 1] Min. number of reads = 456802] Max. number of reads = 4514903] Total number of reads = 60921854] Average number of reads = 121843.75] Median number of reads = 1146517] Sparsity = 0.3219042968756] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 45680"
##
## [[2]]
## [1] "2] Max. number of reads = 451490"
##
## [[3]]
## [1] "3] Total number of reads = 6092185"
##
## [[4]]
## [1] "4] Average number of reads = 121843.7"
##
## [[5]]
## [1] "5] Median number of reads = 114651"
##
## [[6]]
## [1] "7] Sparsity = 0.321904296875"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."Create Ampvis Objects
Create full, core and Site ampvis objects
library(tibble)
library(ampvis2)
# Ampvis object function
create_ampvis_object <- function(pseq) {
# create and extract otu table
otu_table_ampvis <- data.frame(OTU = rownames(phyloseq::otu_table(pseq)@.Data),
phyloseq::otu_table(pseq)@.Data,
phyloseq::tax_table(pseq)@.Data,
check.names = FALSE)
# Metadata
meta_data_ampvis <- data.frame(phyloseq::sample_data(pseq),
check.names = FALSE)
# change index by SampleID
meta_data_ampvis <- meta_data_ampvis %>% rownames_to_column(var = "SampleID")
# ampvis object
av2 <- amp_load(otu_table_ampvis, meta_data_ampvis)
return(av2)
}# phyloseq objects
phyloseq_objects <-list(ps_s1 = ps_s1, ps_s2 = ps_s2, ps_s3 = ps_s3,
ps_s4 = ps_s4, ps_s5 = ps_s5, ps_s6 = ps_s6,
ps_s7 = ps_s7, ps=ps, full_ps_core_counts = full_ps_core_counts)
# Apply function
ampvis_objects <- lapply(phyloseq_objects, create_ampvis_object)
ampvis_objects## $ps_s1
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 6310 1429162 148635 218877 171227
## Avg#Reads
## 178645.25
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 6310(100%) 6291(99.7%) 6261(99.22%) 6055(95.96%) 6044(95.78%) 5946(94.23%)
## Species
## 4349(68.92%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s2
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 5831 1379847 132348 220547 173469.5
## Avg#Reads
## 172480.88
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 5831(100%) 5812(99.67%) 5782(99.16%) 5581(95.71%) 5574(95.59%) 5481(94%)
## Species
## 4018(68.91%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s3
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 4781 1686895 101979 674846 141534
## Avg#Reads
## 210861.88
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 4781(100%) 4764(99.64%) 4735(99.04%) 4557(95.31%) 4552(95.21%) 4477(93.64%)
## Species
## 3322(69.48%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s4
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 4 5408 715609 127631 206504 190737
## Avg#Reads
## 178902.25
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 5408(100%) 5393(99.72%) 5366(99.22%) 5185(95.88%) 5176(95.71%) 5091(94.14%)
## Species
## 3694(68.31%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s5
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 6205 1462992 97081 298225 167958
## Avg#Reads
## 182874
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 6205(100%) 6188(99.73%) 6159(99.26%) 5952(95.92%) 5943(95.78%) 5838(94.09%)
## Species
## 4243(68.38%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s6
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 6 4113 1044090 154919 201962 172266.5
## Avg#Reads
## 174015
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 4113(100%) 4101(99.71%) 4084(99.29%) 3938(95.75%) 3927(95.48%) 3830(93.12%)
## Species
## 2702(65.69%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps_s7
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 8 4860 1287903 66413 215740 174806.5
## Avg#Reads
## 160987.88
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 4860(100%) 4845(99.69%) 4826(99.3%) 4663(95.95%) 4651(95.7%) 4552(93.66%)
## Species
## 3230(66.46%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $ps
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 50 7097 9006498 66413 674846 171486.5
## Avg#Reads
## 180129.96
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 7097(100%) 7077(99.72%) 7043(99.24%) 6811(95.97%) 6795(95.74%) 6672(94.01%)
## Species
## 4818(67.89%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.
##
## $full_ps_core_counts
## ampvis2 object with 3 elements.
## Summary of OTU table:
## Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads
## 50 2048 6092185 45680 451490 114651
## Avg#Reads
## 121843.7
##
## Assigned taxonomy:
## Kingdom Phylum Class Order Family Genus
## 2048(100%) 2043(99.76%) 2031(99.17%) 1960(95.7%) 1957(95.56%) 1913(93.41%)
## Species
## 1336(65.23%)
##
## Metadata variables: 30
## SampleID, Sample, ID, Depth, Depth_cm, Site, Sample_type, Date, Location, Temp...C..int, Salnidad..ppm.int, Salinidad..ppt.int, pH.int, pHsup, Redox.medido.en.campo1, Redox.mV..Eh.int, Temperatura...C.sup, Salnidad..ppm.sup, Salinidad..ppt.sup, Redox.medido.en.campo, Redox.mV..Eh.sup, S.2..mg.l..campo, ml.de.muestra1, Factor.de.dilucion1, S.2..mg.l., ID_M, SO.....mg.l..campo, ml.de.muestra2, Factor.de.dilucion2, SO.....mg.l.#subset amp obj
list(ps_s1 = ps_s1, ps_s2 = ps_s2, ps_s3 = ps_s3,
ps_s4 = ps_s4, ps_s5 = ps_s5, ps_s6 = ps_s6,
ps_s7 = ps_s7, ps=ps)## $ps_s1
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6310 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6310 taxa by 7 taxonomic ranks ]
##
## $ps_s2
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5831 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5831 taxa by 7 taxonomic ranks ]
##
## $ps_s3
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4781 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4781 taxa by 7 taxonomic ranks ]
##
## $ps_s4
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 5408 taxa and 4 samples ]
## sample_data() Sample Data: [ 4 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 5408 taxa by 7 taxonomic ranks ]
##
## $ps_s5
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 6205 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 6205 taxa by 7 taxonomic ranks ]
##
## $ps_s6
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4113 taxa and 6 samples ]
## sample_data() Sample Data: [ 6 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4113 taxa by 7 taxonomic ranks ]
##
## $ps_s7
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 4860 taxa and 8 samples ]
## sample_data() Sample Data: [ 8 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 4860 taxa by 7 taxonomic ranks ]
##
## $ps
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7097 taxa and 50 samples ]
## sample_data() Sample Data: [ 50 samples by 29 sample variables ]
## tax_table() Taxonomy Table: [ 7097 taxa by 7 taxonomic ranks ]av2_full <- ampvis_objects$ps
av2_full_core <- ampvis_objects$full_ps_core_counts
av2_s1 <- ampvis_objects$ps_s1
av2_s2 <-ampvis_objects$ps_s2
av2_s3 <- ampvis_objects$ps_s3
av2_s4 <- ampvis_objects$ps_s4
av2_s5 <- ampvis_objects$ps_s5
av2_s6 <- ampvis_objects$ps_s6
av2_s7 <- ampvis_objects$ps_s7Abundances
Relative abundances
All
#Factor Depth_cm
av2_full$metadata$Depth <- factor(av2_full$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# Full all samples
ampv_heatmap_abundances_phylum_full <- amp_heatmap(av2_full,
group_by = "Depth",
facet_by = "Site",
plot_values = TRUE,
tax_show = 25,
tax_aggregate = "Phylum",
tax_add = "Kingdom",
plot_colorscale = "log10",
plot_legendbreaks = c(1, 5, 5),
color_vector = c("white","deepskyblue3","magenta3"))
ampv_heatmap_abundances_phylum_full
Core
#Factor Depth_cm
av2_full_core$metadata$Depth <- factor(av2_full_core$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# Full all samples
ampv_heatmap_abundances_phylum_core <- amp_heatmap(av2_full_core,
group_by = "Depth",
facet_by = "Site",
plot_values = TRUE,
tax_show = 25,
tax_aggregate = "Phylum",
tax_add = "Kingdom",
plot_colorscale = "log10",
plot_legendbreaks = c(1, 5, 5),
color_vector = c("white","deepskyblue3","magenta3"))
ampv_heatmap_abundances_phylum_core
Combine plots
library(cowplot)
#Combine venn plot
title_abund_plot <- ggdraw() + draw_label("Relative Abundance", fontface = 'bold', x = 0.5, hjust = 0.5, size = 15)
ampv_all_g_plots <- plot_grid(title_abund_plot, ampv_heatmap_abundances_phylum_full, ampv_heatmap_abundances_phylum_core, labels = c("", "A", "B"), ncol = 1, rel_heights = c(0.15, 1, 1))
ampv_all_g_plots
Rank abundance
All
library(ampvis2)
library(ggplot2)
#Factor Depth_cm
av2_full$metadata$Depth <- factor(av2_full$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_full <- amp_rankabundance(av2_full, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
rank_ab_full
#ggsave("results/plots/02.diversity/02.rank_abundance.png", rank_ab , width = 5, height = 4)Site
# Site1
#Factor Depth_cm
av2_s1$metadata$Depth <- factor(av2_s1$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s1 <- amp_rankabundance(av2_s1, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site2
#Factor Depth_cm
av2_s2$metadata$Depth <- factor(av2_s2$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s2 <- amp_rankabundance(av2_s2, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site3
#Factor Depth_cm
av2_s3$metadata$Depth <- factor(av2_s3$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s3 <- amp_rankabundance(av2_s3, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site4
#Factor Depth_cm
av2_s4$metadata$Depth <- factor(av2_s4$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s4 <- amp_rankabundance(av2_s4, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site5
#Factor Depth_cm
av2_s5$metadata$Depth <- factor(av2_s5$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s5 <- amp_rankabundance(av2_s5, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site6
#Factor Depth_cm
av2_s6$metadata$Depth <- factor(av2_s6$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s6 <- amp_rankabundance(av2_s6, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
# Site7
#Factor Depth_cm
av2_s7$metadata$Depth <- factor(av2_s7$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_s7 <- amp_rankabundance(av2_s7, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()Core
library(ampvis2)
library(ggplot2)
#Factor Depth_cm
av2_full_core$metadata$Depth <- factor(av2_full_core$metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# rank abundance
rank_ab_full_core <- amp_rankabundance(av2_full_core, log10_x = TRUE, group_by = "Depth") +
scale_color_lancet()
rank_ab_full_core
Combine
library(cowplot)
#Combine euler plot
title_rank_plot <- ggdraw() + draw_label("Relative diversity and species evenness", fontface = 'bold', x = 0.5, hjust = 0.5, size = 15)
rank_ab_full_core_s1 <- plot_grid(rank_ab_full, rank_ab_full_core, rank_ab_s1, labels = c("F", "C", "S1"), ncol = 3)
rank_ab_s2tos4 <- plot_grid(rank_ab_s2,rank_ab_s3, rank_ab_s4, labels = c("S2","S3", "S4"), ncol = 3)## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Infrank_ab_s5tos7 <- plot_grid(rank_ab_s5, rank_ab_s6, rank_ab_s7, labels = c("S5", "S6", "S7"), ncol = 3)
rankabund_all_plots <- plot_grid(title_rank_plot, rank_ab_full_core_s1, rank_ab_s2tos4, rank_ab_s5tos7, ncol = 1, rel_heights = c(0.14, 1, 1, 1))
rankabund_all_plots
Alpha diversity
Full
Get Hill numbers
library(hilldiv)## Registered S3 methods overwritten by 'FSA':
## method from
## confint.boot car
## hist.boot carlibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ lubridate 1.9.3 ✔ stringr 1.5.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ microbiome::alpha() masks ggplot2::alpha()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ lubridate::stamp() masks cowplot::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(kableExtra)##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows# Get hill numbers
q0 <- hill_div(otu_table, qvalue = 0)
q1 <- hill_div(otu_table, qvalue = 1)
q2 <- hill_div(otu_table, qvalue = 2)
# Merge metadata with Hill numbers
q012_all <- cbind(q0, q1, q2) %>% as.data.frame() %>% rownames_to_column(var = "SampleID")
metadata_with_hill <- q012_all %>%
inner_join(metadata, by = c("SampleID"="SampleID"))
#save table
write.table(metadata_with_hill, "../results/05.diversity/Metadata_with_hill.tsv", quote = FALSE, sep = "\t", row.names = TRUE, col.names=TRUE)
# show
metadata_with_hill %>% head() %>% kable()| SampleID | q0 | q1 | q2 | Sample | ID | Depth | Depth_cm | Site | Sample_type | Date | Location | Temp (°C) int | Salnidad (ppm)int | Salinidad (ppt)int | pH int | pHsup | Redox medido en campo1 | Redox mV (Eh)int | Temperatura (°C)sup | Salnidad (ppm)sup | Salinidad (ppt)sup | Redox medido en campo | Redox mV (Eh)sup | S-2 (mg/l) campo | ml de muestra1 | Factor de dilucion1 | S-2 (mg/l) | ID_M | SO₄²- (mg/l) campo | ml de muestra2 | Factor de dilucion2 | SO₄ ² (mg/l) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| S1A015 | 2267 | 581.0228 | 172.50572 | S1A015 | S1A | 1 | 0-15 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1A1530 | 3082 | 924.7261 | 341.75340 | S1A1530 | S1A | 2 | 16-30 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1A3050 | 3234 | 926.0288 | 303.32359 | S1A3050 | S1A | 3 | 31-45 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1A5075 | 3506 | 879.5831 | 228.44914 | S1A5075 | S1A | 4 | 50-75 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 31.4 | 817.5 | 0.82 | 6.55 | 7.34 | -114 | 121 | 33.0 | 734.7 | 0.73 | -108 | 127 | 7.73 | 25 | 1 | 7.73 | S1a M | 53.7 | 2 | 10 | 537 |
| S1B07 | 3011 | 364.4166 | 62.98918 | S1B07 | S1B | 1 | 0-15 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 30.9 | 794.0 | 0.79 | 6.69 | 7.79 | -234 | 1 | 33.2 | 732.8 | 0.73 | -136 | 99 | 7.09 | 25 | 1 | 7.09 | S1b M | 55.0 | 2 | 10 | 550 |
| S1B2230 | 3973 | 1127.1162 | 320.30757 | S1B2230 | S1B | 2 | 16-30 | Site_1 | Sediment | 28/09/2023 | Dique Miguelito | 30.9 | 794.0 | 0.79 | 6.69 | 7.79 | -234 | 1 | 33.2 | 732.8 | 0.73 | -136 | 99 | 7.09 | 25 | 1 | 7.09 | S1b M | 55.0 | 2 | 10 | 550 |
Explore Hill distribution
Get means
# Reorder q values
meta_qs <- metadata_with_hill %>%
pivot_longer(cols = q0:q2, names_to = "q", values_to = "value") %>%
filter(q %in% c("q0", "q1", "q2")) %>%
mutate(
qs = case_when(
q == "q0" ~ "q=0",
q == "q1" ~ "q=1",
q == "q2" ~ "q=2"
))
#Get means of Hill numbers
means <- meta_qs %>% group_by(Depth_cm, Site, qs) %>%
summarise(means = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
.groups = 'drop')
#save table
write.table(means, "../results/05.diversity/Hill_means_sd.tsv", quote = FALSE, sep = "\t", row.names = TRUE, col.names=TRUE)
print(means)## # A tibble: 81 × 5
## Depth_cm Site qs means sd
## <chr> <chr> <chr> <dbl> <dbl>
## 1 0-15 Site_1 q=0 2639 526.
## 2 0-15 Site_1 q=1 473. 153.
## 3 0-15 Site_1 q=2 118. 77.4
## 4 0-15 Site_2 q=0 3142. 771.
## 5 0-15 Site_2 q=1 941. 339.
## 6 0-15 Site_2 q=2 288. 116.
## 7 0-15 Site_3 q=0 2771 781.
## 8 0-15 Site_3 q=1 612. 429.
## 9 0-15 Site_3 q=2 150. 123.
## 10 0-15 Site_4 q=0 3009 NA
## # ℹ 71 more rowsPlot Hill distribution
library(ggpubr) ##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legend## The following object is masked from 'package:qiime2R':
##
## mean_sd# factor to reorder plot
metadata_with_hill$Site <- factor(metadata_with_hill$Site)
metadata_with_hill$Depth <- factor(metadata_with_hill$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
#plot
hill_barplot_explore <- metadata_with_hill %>%
pivot_longer(cols = q0:q2, names_to = "q", values_to = "value") %>%
filter(q %in% c("q0", "q1", "q2")) %>%
mutate(
qs = case_when(
q == "q0" ~ "q=0",
q == "q1" ~ "q=1",
q == "q2" ~ "q=2"
)) %>%
ggbarplot(
x = "Depth",
y = "value",
add = "mean_se",
facet.by = c("qs", "Site"),
fill = "Depth"
) +
scale_fill_lancet() +
#geom_jitter(size = 0.8, position = position_jitter(width = 0.1)) +
facet_grid(rows = vars(qs), cols = vars(Site), scales = "free_y") +
theme(
strip.background = element_blank(), #element_rect(fill = "")
strip.text.x = element_text(face= "bold", size = 13, margin = margin(0.5, 0, 0.5, 0)),
strip.text.y = element_text(face= "bold", size = 14, angle = 0, margin = margin(0, 0.5, 0, 0.5)),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line = element_line(colour = "black"),
#panel.border = element_blank(),
panel.spacing.x = unit(0.5, "lines"),
panel.background = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_text(colour = "black", size = 11),
legend.position = "bottom"
) +
labs(
fill = "Depth",
y = "",
x = "",
title = "Biodiversity across Sites and Depths"
) +
theme(legend.text = element_text(size = 12))## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs produced
## Warning in stats::qt(ci/2 + 0.5, data_sum$length - 1): NaNs producedhill_barplot_explore
#ggsave("../results/plots/Hill_diversity_across_Sites_and_Depths_fltr.png", hill_barplot_explore_fltr, width = 11, height = 11)Depth vs Hill Correlation
Alpha diversity depth correlation to order q=0
#Get effective reads
Effective_reads <- colSums(otu_table[-1, ])
metadata_with_hill$Effective_reads <- Effective_reads#install.packages("ggpubr")
library(ggpubr)
q0_vs_depth <- ggscatter(metadata_with_hill,
x = "Effective_reads",
y = "q0",
xlab= "Sequencing depth (number of reads)",
ylab = "q=0",
#ylab="Alpha diversity q=0 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#114e9d", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")Alpha diversity depth correlation to order q=1
q1_vs_depth <- ggscatter(metadata_with_hill,
x = "Effective_reads",
y = "q1",
xlab= "Sequencing depth (number of reads)",
ylab = "q=1",
#ylab="Alpha diversity q=1 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#114e9d", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")Alpha diversity depth correlation to order q=2
q2_vs_depth <- ggscatter(metadata_with_hill,
x = "Effective_reads",
y = "q2",
xlab= "Sequencing depth (number of reads)",
ylab = "q=2",
#ylab="Alpha diversity q=2 (effective number of total ASVs)",
add = "reg.line", # Add regression line
add.params = list(color = "#114e9d", fill = "lightgray"), # Customize reg. line
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")
)+
theme(legend.title = element_blank(), legend.position = "none")library(ggplot2)
library(cowplot)
#save plots
#Combine plot
title_corr_plot <- ggdraw() + draw_label("Alpha diversity depth correlation with fltr samples")
correlation_plot_q012 <- plot_grid(title_corr_plot, q0_vs_depth, q1_vs_depth, q2_vs_depth, labels = c(" ","A", "B", "C"), ncol = 1, rel_heights = c(0.15, 1, 1, 1))
correlation_plot_q012
#save plot
#ggsave("../results/plots/alpha_depth_correlation_samples.png", correlation_plot_q012_, width = 8, height = 8)Check Hill numbers Normality
Shapiro test to Hill numbers
library(ggplot2)
library(cowplot)
#Shapiro test
shapiro_q0 <- shapiro.test(metadata_with_hill$q0)
shapiro_q0_pvalue <- round(shapiro_q0$p.value, 5)
shapiro_q1 <- shapiro.test(metadata_with_hill$q1)
shapiro_q1_pvalue <- round(shapiro_q1$p.value, 5)
shapiro_q2 <- shapiro.test(metadata_with_hill$q2)
shapiro_q2_pvalue <- round(shapiro_q2$p.value, 5)
#Histograms
histplot_q0 <- ggplot(metadata_with_hill, aes(x = q0, xlab="q=0")) +
geom_histogram(fill = "lightblue", bins = 15) +
ggtitle(paste("Shapiro, p-value:", shapiro_q0_pvalue)) +
theme_bw() + xlab("q=0") + ylab("Frequency")
histplot_q1 <- ggplot(metadata_with_hill, aes(x = q1)) +
geom_histogram(fill = "lightblue", bins = 15) +
ggtitle(paste("Shapiro, p-value:", shapiro_q1_pvalue)) +
theme_bw() + xlab("q=1") + ylab("Frequency")
histplot_q2 <- ggplot(metadata_with_hill, aes(x = q2)) +
geom_histogram(fill = "lightblue", bins = 15) +
ggtitle(paste("Shapiro, p-value:", shapiro_q2_pvalue)) +
theme_bw() + xlab("q=2") + ylab("Frequency")
#Combine plot
title_plot <- ggdraw() + draw_label("Histogram of Hill diversity") #, fontface = 'bold', x = 0.5, hjust = 0.5)
histplot_q012 <- plot_grid(title_plot, histplot_q0, histplot_q1, histplot_q2, labels = c(" ","A", "B", "C"), ncol = 1, rel_heights = c(0.15, 1, 1, 1))
histplot_q012
#save plot
#ggsave("../results/plots/hill_normality_test_histogram_samples.png", histplot_q012)A pesar de que tienen una distribución normal, dado que los datos del proyecto tienen medidas repetidas y desbalanceo, lo adecuado es hacer comparaciones basadas en modelos. En este caso modelos lineales mixtos
Models
Linear Mixed Models
Estos modelos son útiles cuando las observaciones tienen estructuras de dependencia, como datos longitudinales o datos agrupados/anidados. Los modelos lineales mixtos incluyen tanto efectos fijos como efectos aleatorios, donde los efectos fijos son comunes a todas las unidades/individuos, y los efectos aleatorios varÃan por grupo o unidad.
lme
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapselibrary(pgirmess)## Registered S3 method overwritten by 'spdep':
## method from
## plot.mst ape# lme
# q0
q0_lme <- lme(q0 ~ log(Effective_reads)+Depth_cm+Site, random=~1 |ID, data = metadata_with_hill)
q0_lme_perm <- PermTest(q0_lme, B = 100)
# q1
q1_lme <- lme(q1 ~ log(Effective_reads)+Depth_cm+Location, random=~1 |ID, data = metadata_with_hill)
q1_lme_perm <- PermTest(q1_lme, B = 100)
# q2
q2_lme <- lme(q2 ~ log(Effective_reads)+Depth_cm+Location, random=~1 |ID, data = metadata_with_hill)
q2_lme_perm <- PermTest(q2_lme, B = 100)Effect
# extract p-values fun
extract_pvalues <- function(perm_test_result) {
perm_test_result$resultats$p.value
}
# extract
lme_pvalues_q0 <- extract_pvalues(q0_lme_perm)
lme_pvalues_q1 <- extract_pvalues(q1_lme_perm)
lme_pvalues_q2 <- extract_pvalues(q2_lme_perm)
# Dataframe p-values
library(dplyr)
lme_pvalues_table <- data.frame(
Hill = c("q0", "q1", "q2"),
`Intercept` = c(lme_pvalues_q0[1], lme_pvalues_q1[1], lme_pvalues_q2[1]),
`log(Depth)` = c(lme_pvalues_q0[2], lme_pvalues_q1[2], lme_pvalues_q2[2]),
`Depth` = c(lme_pvalues_q0[3], lme_pvalues_q1[3], lme_pvalues_q2[3]),
`Site` = c(lme_pvalues_q0[4], lme_pvalues_q1[4], lme_pvalues_q2[4])
)
# table
knitr::kable(t(lme_pvalues_table), caption = "P-values from Permuted LME")| Hill | q0 | q1 | q2 |
| Intercept | 0.17 | 0.63 | 0.61 |
| log.Depth. | 0.00 | 0.95 | 0.81 |
| Depth | 0.00 | 0.14 | 0.08 |
| Site | 0.20 | 0.81 | 0.99 |
# save
write.csv(t(lme_pvalues_table), "../results/05.diversity/P-values_from_Permuted_LME.csv", row.names = FALSE)q0_lme_perm##
## Monte-Carlo test
##
## Call:
## PermTest.lme(obj = q0_lme, B = 100)
##
## Based on 100 replicates
## Simulated p-value:
## p.value
## (Intercept) 0.17
## log(Effective_reads) 0.00
## Depth_cm 0.00
## Site 0.20Check with easystats
q=0
library(easystats)## # Attaching packages: easystats 0.7.1 (red = needs update)
## ✔ bayestestR 0.13.2 ✔ correlation 0.8.4
## ✔ datawizard 0.10.0 ✔ effectsize 0.8.7
## ✔ insight 0.19.10 ✔ modelbased 0.8.7
## ✔ performance 0.11.0 ✔ parameters 0.21.6
## ✔ report 0.5.8 ✖ see 0.8.3
##
## Restart the R-Session and update packages with `easystats::easystats_update()`.check_q0_lme_1 <- check_model(q0_lme)## Converting missing values (`NA`) into regular values currently not
## possible for variables of class `NULL`.print(check_q0_lme_1)
q=1
check_q1_lme_1 <- check_model(q1_lme)## Converting missing values (`NA`) into regular values currently not
## possible for variables of class `NULL`.print(check_q1_lme_1)
q=2
check_q2_lme_1 <- check_model(q2_lme)## Converting missing values (`NA`) into regular values currently not
## possible for variables of class `NULL`.check_q2_lme_1
Generalized Linear Mixed Model
glm quasipoisson
Beta diversity Diversity
Full
clr manual
Prepare data
aldex_clr_data <- t(getMonteCarloSample(aldex_clr,1))
pca <- prcomp(aldex_clr_data)
meta=metadata_with_hill %>% dplyr::rename(SampleID="SampleID")PCA
library(ggplot2)
library(dplyr)
library(tidyr)
# Reorder
meta$Depth <- factor(meta$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
# plot
pca_plot <- ggplot(data = data.frame(pca$x) %>%
rownames_to_column(var = "SampleID") %>%
left_join(meta, by = "SampleID")) +
geom_segment(data = data.frame(pca$rotation) %>%
rownames_to_column(var = "Feature.ID") %>%
mutate(a = sqrt(PC1^2 + PC2^2)) %>%
top_n(8, a) %>%
mutate(PC1 = PC1*50, PC2 = PC2*50),
aes(x = 0, xend = PC1, y = 0, yend = PC2)) +
geom_point(aes(x = PC1, y = PC2, color = Depth, shape = Site), size = 3) +
geom_vline(xintercept = 0, linetype = 2) +
geom_hline(yintercept = 0, linetype = 2) +
scale_color_lancet() +
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
scale_fill_lancet() +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical") +
labs(y = "PC2", x = "PC1", color = "Depth", shape = "Site") +
theme_bw()
pca_plot
bray
# Bray NMDS bray
library(vegan)## Loading required package: permute##
## Attaching package: 'permute'## The following object is masked from 'package:devtools':
##
## check## Loading required package: lattice## This is vegan 2.6-4##
## Attaching package: 'vegan'## The following object is masked from 'package:microbiome':
##
## diversitybray=vegdist(t(otu_table), method = "bray")
nmds_source_bray = vegan::metaMDS(bray, trymax = 20, k = 2) ## Run 0 stress 0.1712067
## Run 1 stress 0.1756919
## Run 2 stress 0.1712066
## ... New best solution
## ... Procrustes: rmse 0.0001310485 max resid 0.0008477333
## ... Similar to previous best
## Run 3 stress 0.1742446
## Run 4 stress 0.1850122
## Run 5 stress 0.1885669
## Run 6 stress 0.1712066
## ... Procrustes: rmse 0.00003115916 max resid 0.0001995856
## ... Similar to previous best
## Run 7 stress 0.1712066
## ... Procrustes: rmse 0.00001631054 max resid 0.0001018519
## ... Similar to previous best
## Run 8 stress 0.1946932
## Run 9 stress 0.1742446
## Run 10 stress 0.1756918
## Run 11 stress 0.1756918
## Run 12 stress 0.1712068
## ... Procrustes: rmse 0.0002026832 max resid 0.001330524
## ... Similar to previous best
## Run 13 stress 0.1874831
## Run 14 stress 0.1742447
## Run 15 stress 0.1850122
## Run 16 stress 0.1900205
## Run 17 stress 0.1757299
## Run 18 stress 0.1741713
## Run 19 stress 0.1777425
## Run 20 stress 0.1712067
## ... Procrustes: rmse 0.000111613 max resid 0.0007341006
## ... Similar to previous best
## *** Best solution repeated 5 timesvar_stress_nmds_bray <- round(nmds_source_bray$stress, 5)
var_stress_nmds_bray## [1] 0.17121scores_source_bray= nmds_source_bray %>% vegan::scores()
metadata$Depth <- factor(metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
nmds_plot_bray<- ggplot() +
geom_point(data=data.frame(scores_source_bray) %>%
rownames_to_column(var = "SampleID")%>%
left_join(metadata, by = "SampleID"),
aes(x=NMDS1, y=NMDS2, color = Depth, shape = Site), size=5) +
theme_linedraw()+
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw() +
labs(title="NMDS, Bray-Curtis distances") +
ylab("NMDS2")+xlab("NMDS1")
nmds_plot_bray <- nmds_plot_bray +
annotate("text", x = Inf, y = -Inf, label = paste("Stress:", var_stress_nmds_bray),
hjust = 3.1, vjust = -1.9, size = 3)
nmds_plot_bray
aitchison
# NMDS aichitson
library(vegan)
pseudocount <- 1
otu_table_pseudo <- otu_table + pseudocount
aitch=vegdist(t(otu_table_pseudo), method = "aitchison")
nmds_source_aitch = vegan::metaMDS(aitch, trymax = 20, k = 2) ## Run 0 stress 0.1034177
## Run 1 stress 0.1034177
## ... Procrustes: rmse 0.00000122641 max resid 0.000003746513
## ... Similar to previous best
## Run 2 stress 0.1542678
## Run 3 stress 0.1519305
## Run 4 stress 0.1034177
## ... Procrustes: rmse 0.000005404249 max resid 0.00002283729
## ... Similar to previous best
## Run 5 stress 0.1034177
## ... Procrustes: rmse 0.000001584979 max resid 0.000005991313
## ... Similar to previous best
## Run 6 stress 0.1519474
## Run 7 stress 0.1056559
## Run 8 stress 0.1056635
## Run 9 stress 0.1080944
## Run 10 stress 0.1034177
## ... Procrustes: rmse 0.000002979132 max resid 0.00001396478
## ... Similar to previous best
## Run 11 stress 0.1034177
## ... New best solution
## ... Procrustes: rmse 0.000001200118 max resid 0.000003370918
## ... Similar to previous best
## Run 12 stress 0.1034177
## ... Procrustes: rmse 0.000001639774 max resid 0.000004650393
## ... Similar to previous best
## Run 13 stress 0.1080944
## Run 14 stress 0.1034177
## ... Procrustes: rmse 0.000001682944 max resid 0.000006160126
## ... Similar to previous best
## Run 15 stress 0.1056635
## Run 16 stress 0.1080944
## Run 17 stress 0.1034177
## ... New best solution
## ... Procrustes: rmse 0.00000116747 max resid 0.000003465607
## ... Similar to previous best
## Run 18 stress 0.1034177
## ... Procrustes: rmse 0.000001206626 max resid 0.000002998858
## ... Similar to previous best
## Run 19 stress 0.1080944
## Run 20 stress 0.1034177
## ... Procrustes: rmse 0.000001326051 max resid 0.000004422972
## ... Similar to previous best
## *** Best solution repeated 3 timesvar_stress_nmds_aitch <- round(nmds_source_aitch$stress, 5)
var_stress_nmds_aitch## [1] 0.10342scores_source_aitch= nmds_source_aitch %>% vegan::scores()
metadata$Depth <- factor(metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
nmds_plot_aitch<- ggplot() +
geom_point(data=data.frame(scores_source_aitch) %>%
rownames_to_column(var = "SampleID")%>%
left_join(metadata, by = "SampleID"),
aes(x=NMDS1, y=NMDS2, color = Depth, shape = Site), size=5) +
# theme_linedraw()+
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw() +
labs(title="NMDS, Aitchison distances") +
ylab("NMDS2")+xlab("NMDS1")
nmds_plot_aitch <- nmds_plot_aitch +
annotate("text", x = Inf, y = -Inf, label = paste("Stress:", var_stress_nmds_aitch),
hjust = 3.1, vjust = -1.9, size = 3)
nmds_plot_aitch
library(vegan)
bray_dist =vegdist(t(otu_table), method = "bray")
pcoas=wcmdscale(bray_dist, eig = T)
metadata$Depth <- factor(metadata$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
pcoa_plot<- ggplot() +
geom_point(data=data.frame(pcoas$points) %>% #individuals
rownames_to_column(var = "SampleID")%>%
left_join(metadata, by = "SampleID"),
aes(x=Dim1, y=Dim2,color = Depth, shape = Site ), size=4) +
geom_vline(xintercept = 0, linetype = 2) + #lines-cross
geom_hline(yintercept = 0, linetype = 2) +
theme_linedraw()+
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw()+
labs(title="PCoA, Bray-Curtis distances")+
ylab("PCoA2")+xlab("PCoA1")
pcoa_plot
#Combine plot
nmds_bray_aitch <- plot_grid(nmds_plot_bray, nmds_plot_aitch, labels = c("A", "B"), ncol = 2)
nmds_bray_aitch
Permanova
meta_just = data.frame(t(otu_table_pseudo), check.names = F) %>% rownames_to_column(var = "SampleID") %>% inner_join(metadata_with_hill)## Joining with `by = join_by(SampleID)`permanovas = adonis2(t(otu_table_pseudo) ~ log(Effective_reads)+Depth_cm*Site,
data=meta_just, method = "aitchison",
permutations = 999, strata = meta_just$ID)
permanovasCore
Aitchison
Prepare data
Explore
microbiome::summarize_phyloseq(full_ps_core_counts)## Compositional = NO2## 1] Min. number of reads = 456802] Max. number of reads = 4514903] Total number of reads = 60921854] Average number of reads = 121843.75] Median number of reads = 1146517] Sparsity = 0.3219042968756] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons
## (i.e. exactly one read detected across all samples)010] Number of sample variables are: 29SampleIDDepthDepth_cmSiteSample_typeDateLocationTemp...C..intSalnidad..ppm.intSalinidad..ppt.intpH.intpHsupRedox.medido.en.campo1Redox.mV..Eh.intTemperatura...C.supSalnidad..ppm.supSalinidad..ppt.supRedox.medido.en.campoRedox.mV..Eh.supS.2..mg.l..campoml.de.muestra1Factor.de.dilucion1S.2..mg.l.ID_MSO.....mg.l..campoml.de.muestra2Factor.de.dilucion2SO.....mg.l.2## [[1]]
## [1] "1] Min. number of reads = 45680"
##
## [[2]]
## [1] "2] Max. number of reads = 451490"
##
## [[3]]
## [1] "3] Total number of reads = 6092185"
##
## [[4]]
## [1] "4] Average number of reads = 121843.7"
##
## [[5]]
## [1] "5] Median number of reads = 114651"
##
## [[6]]
## [1] "7] Sparsity = 0.321904296875"
##
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
##
## [[8]]
## [1] "8] Number of singletons = 0"
##
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n (i.e. exactly one read detected across all samples)0"
##
## [[10]]
## [1] "10] Number of sample variables are: 29"
##
## [[11]]
## [1] "Sample" "ID" "Depth"
## [4] "Depth_cm" "Site" "Sample_type"
## [7] "Date" "Location" "Temp...C..int"
## [10] "Salnidad..ppm.int" "Salinidad..ppt.int" "pH.int"
## [13] "pHsup" "Redox.medido.en.campo1" "Redox.mV..Eh.int"
## [16] "Temperatura...C.sup" "Salnidad..ppm.sup" "Salinidad..ppt.sup"
## [19] "Redox.medido.en.campo" "Redox.mV..Eh.sup" "S.2..mg.l..campo"
## [22] "ml.de.muestra1" "Factor.de.dilucion1" "S.2..mg.l."
## [25] "ID_M" "SO.....mg.l..campo" "ml.de.muestra2"
## [28] "Factor.de.dilucion2" "SO.....mg.l."#Get otu and metadata table of core phyloseq object
# extract otu table
otu_table_full_core <- otu_table(full_ps_core_counts@otu_table)
meta_full_core <- sample_data(full_ps_core_counts)# NMDS aichitson
library(vegan)
pseudocount <- 1
otu_table_full_core_pseudo <- otu_table_full_core + pseudocount
aitch_fc=vegdist(t(otu_table_full_core_pseudo), method = "aitchison")
nmds_source_aitch_fc = vegan::metaMDS(aitch_fc, trymax = 20, k = 2) ## Run 0 stress 0.09967979
## Run 1 stress 0.09967979
## ... Procrustes: rmse 0.000007858841 max resid 0.0000363496
## ... Similar to previous best
## Run 2 stress 0.1509729
## Run 3 stress 0.1023438
## Run 4 stress 0.09967979
## ... New best solution
## ... Procrustes: rmse 0.000004178664 max resid 0.0000199997
## ... Similar to previous best
## Run 5 stress 0.09967979
## ... Procrustes: rmse 0.000005413533 max resid 0.00002386895
## ... Similar to previous best
## Run 6 stress 0.1023717
## Run 7 stress 0.1023438
## Run 8 stress 0.1023438
## Run 9 stress 0.1722353
## Run 10 stress 0.09967979
## ... Procrustes: rmse 0.000001449633 max resid 0.000004770371
## ... Similar to previous best
## Run 11 stress 0.09967979
## ... Procrustes: rmse 0.000004391546 max resid 0.0000213067
## ... Similar to previous best
## Run 12 stress 0.09967979
## ... Procrustes: rmse 0.000002657153 max resid 0.00001148891
## ... Similar to previous best
## Run 13 stress 0.1023438
## Run 14 stress 0.09967979
## ... Procrustes: rmse 0.000001907427 max resid 0.000007322214
## ... Similar to previous best
## Run 15 stress 0.09967979
## ... Procrustes: rmse 0.000006794376 max resid 0.00003191459
## ... Similar to previous best
## Run 16 stress 0.1509729
## Run 17 stress 0.1023717
## Run 18 stress 0.1968676
## Run 19 stress 0.09967979
## ... Procrustes: rmse 0.000003876331 max resid 0.0000177978
## ... Similar to previous best
## Run 20 stress 0.1509889
## *** Best solution repeated 8 timesvar_stress_nmds_aitch_fc <- round(nmds_source_aitch_fc$stress, 5)
var_stress_nmds_aitch_fc## [1] 0.09968scores_source_aitch_fc= nmds_source_aitch_fc %>% vegan::scores()
meta_full_core$Depth <- factor(meta_full_core$Depth_cm, levels = c("0-15", "16-30", "31-45", "50-75"))
nmds_plot_aitch_fc<- ggplot() +
geom_point(data=data.frame(scores_source_aitch_fc) %>%
rownames_to_column(var = "Sample")%>%
left_join(meta_full_core, by = "Sample"),
aes(x=NMDS1, y=NMDS2, color = Depth, shape = Site), size=4) +
scale_fill_lancet()+
scale_color_lancet()+
scale_shape_manual(values = c(15, 16, 17, 18, 23, 24, 25)) +
theme(axis.text = element_text(colour = "black", size = 12),
axis.title = element_text(colour = "black", size = 12),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right",
legend.box = "vertical")+ theme_bw() +
labs(title="NMDS, Aitchison distances Core") +
ylab("NMDS2")+xlab("NMDS1")
nmds_plot_aitch_fc <- nmds_plot_aitch_fc +
annotate("text", x = Inf, y = -Inf, label = paste("Stress:", var_stress_nmds_aitch_fc),
hjust = 3.1, vjust = -1.9, size = 3)
nmds_plot_aitch_fc
#Combine plot
nmds_aitch_full_core <- plot_grid(nmds_plot_aitch, nmds_plot_aitch_fc, labels = c("A", "B"), ncol = 2)
nmds_aitch_full_core
Permanova
pseudocount <- 1
otu_table_full_core_pseudo <- otu_table_full_core + pseudocount
#Get effective reads
Effective_reads_core <- colSums(otu_table_full_core[-1, ])
meta_full_core$Effective_reads <- Effective_reads_coremeta_core_pseudo = data.frame(t(otu_table_full_core_pseudo), check.names = F) %>%
rownames_to_column(var = "Sample") %>% inner_join(meta_full_core) ## Joining with `by = join_by(Sample)`perm_core = adonis2(t(otu_table_full_core_pseudo) ~ log(Effective_reads)+Depth_cm*Site,
data=meta_core_pseudo, method = "aitchison",
permutations = 999, strata = meta_core_pseudo$ID)
perm_coreCore Microbiome
Microbiome
Based here
Calculate core between Depths
# count number of samples in each group
table(meta(ps)$Depth_cm, useNA = "always")##
## 0-15 16-30 31-45 50-75 <NA>
## 13 13 13 11 0Full Core Depth
#List of Depths
Depths <- unique(as.character(meta(ps_full_rel_f)$Depth_cm))
print(Depths)## [1] "0-15" "16-30" "31-45" "50-75"#Write a for loop to go through each of the Depths one by one and combine identified core taxa into a list.
list_core_full <- list()# an empty object to store information
for (Depth_cm in Depths){
#print(paste0("Identifying Core Taxa for ", Depth))
ps_sub_full <- subset_samples(ps_full_rel_f, Depth_cm == Depth_cm )
core_m_full <- core_members(ps_sub_full,
detection = 0, # 100 % samples
prevalence = 0.5) # 50% prevalence
print(paste0("No. of core taxa in ", Depth_cm, " : ", length(core_m_full)))
# add to a list core taxa for each group.
list_core_full[[Depth_cm]] <- core_m_full
#create core depth variable
variable_name_full <- paste0("core_full_", gsub(" ", "_", Depth_cm))
assign(variable_name_full, core_m_full, envir = .GlobalEnv)
#write all cores
file_name_full <- paste0("../results/05.diversity/",variable_name_full,".csv")
write.csv(core_m_full, file = file_name_full, row.names = FALSE)
}## [1] "No. of core taxa in 0-15 : 2048"
## [1] "No. of core taxa in 16-30 : 2048"
## [1] "No. of core taxa in 31-45 : 2048"
## [1] "No. of core taxa in 50-75 : 2048"Plot venn diagram
library(VennDiagram)## Loading required package: grid## Loading required package: futile.logger##
## Attaching package: 'VennDiagram'## The following object is masked from 'package:ggpubr':
##
## rotate# Extract lancet Depth colors
depth_colors <- scale_fill_lancet()$palette(4)
# plot
venn_plot_full <- venn.diagram(
x = list_core_full,
filename = NULL,
#filename = "../results/05.diversity/Venn_full.png",
category.names = c("0-15", "16-30", "31-45", "50-75"),
output = TRUE,
#imagetype = "png",
height = 2000,
width = 2000,
resolution = 300,
compression = "lzw",
#units = "px",
lwd = 2,
lty = 'blank',
fontfamily ="sans",
fill = depth_colors,
cat.fontfamily = "sans",
cat.cex = 1.85,
cat.default.pos = "outer")
# Show
grid.draw(venn_plot_full)
library(UpSetR)##
## Attaching package: 'UpSetR'## The following object is masked from 'package:lattice':
##
## histogram# Convert UpSetR fromList
full_core_data <- fromList(list_core_full)
# UpSet plot
upset_full_depth <- upset(full_core_data, nsets = 4, decreasing = T, order.by = "freq", keep.order = TRUE, main.bar.color = "black", sets.bar.color = depth_colors)
upset_full_depth
Get core identities
Full Core Site
#List of Depths
sites <- unique(as.character(meta(ps_full_rel_f)$Site))
print(sites)## [1] "Site_1" "Site_2" "Site_3" "Site_4" "Site_5" "Site_6" "Site_7"#Write a for loop to go through each of the Depths one by one and combine identified core taxa into a list.
list_core_sites_full <- list()# an empty object to store information
for (site in sites){
#print(paste0("Identifying Core Taxa for ", Depth))
ps_sub_full <- subset_samples(ps_full_rel_f, Site == site)
core_m_full <- core_members(ps_sub_full,
detection = 0, # 100 % samples
prevalence = 0.5) # 50% prevalence
print(paste0("No. of core taxa in ", site, " : ", length(core_m_full)))
# add to a list core taxa for each group.
list_core_sites_full[[site]] <- core_m_full
#create core depth variable
variable_name_full <- paste0("core_full_", gsub(" ", "_", site))
assign(variable_name_full, core_m_full, envir = .GlobalEnv)
#write all cores
file_name_full <- paste0("../results/05.diversity/",variable_name_full,".csv")
write.csv(core_m_full, file = file_name_full, row.names = FALSE)
}## [1] "No. of core taxa in Site_1 : 2403"
## [1] "No. of core taxa in Site_2 : 3356"
## [1] "No. of core taxa in Site_3 : 2192"
## [1] "No. of core taxa in Site_4 : 3238"
## [1] "No. of core taxa in Site_5 : 2900"
## [1] "No. of core taxa in Site_6 : 2158"
## [1] "No. of core taxa in Site_7 : 2809"Upset diagram
library(UpSetR)
sites_colors <- scale_fill_d3()$palette(7)
# Convert UpSetR fromList
input_data <- fromList(list_core_sites_full)
# UpSet plot
upset_full_sites <- upset(input_data, nsets = 7, decreasing = T, order.by = "freq", keep.order = TRUE, main.bar.color = "black", sets.bar.color = sites_colors)
upset_full_sites